"""
数据预处理模块 - 阶段1: 初始化与预处理
实现DNI数据处理、太阳轨迹计算、可行域生成等功能
"""

import numpy as np
import pandas as pd
from datetime import datetime, timedelta
import math

class DataPreprocessor:
    def __init__(self):
        """初始化数据预处理器"""
        # 基础参数设置
        self.latitude = 39.4  # 北京纬度 (度)
        self.longitude = 116.4  # 北京经度 (度) 
        self.timezone = 8  # 东八区
        
        # 场地参数
        self.R_field = 350  # 场地半径 (m)
        self.R_keepout = 100  # 禁区半径 (m)
        
        # 定日镜参数 (问题1给定)
        self.mirror_width = 6  # 定日镜宽度 (m)
        self.mirror_height = 6  # 定日镜高度 (m)
        self.mirror_area = self.mirror_width * self.mirror_height  # 镜面面积 (m²)
        self.install_height = 4  # 安装高度 (m)
        
        # 吸收塔参数
        self.tower_height = 80  # 吸收塔高度 (m)
        self.receiver_height = 8  # 接收器高度 (m)
        self.receiver_width = 7  # 接收器宽度 (m)
        
        print("✅ 数据预处理器初始化完成")
    
    def load_boundary_data(self, filename='附件.xlsx'):
        """加载场地边界数据"""
        try:
            self.boundary_data = pd.read_excel(filename)
            x_coords = self.boundary_data.iloc[:, 0].values
            y_coords = self.boundary_data.iloc[:, 1].values
            
            print(f"✅ 场地边界数据加载成功: {len(x_coords)} 个边界点")
            print(f"   X坐标范围: [{x_coords.min():.1f}, {x_coords.max():.1f}] m")
            print(f"   Y坐标范围: [{y_coords.min():.1f}, {y_coords.max():.1f}] m")
            
            return x_coords, y_coords
        except Exception as e:
            print(f"❌ 场地边界数据加载失败: {e}")
            return None, None
    
    def generate_time_points(self, year=2024):
        """生成全年时间点 (每小时采样)"""
        time_points = []
        start_date = datetime(year, 1, 1, 0, 0, 0)
        
        # 生成全年每小时的时间点
        for day in range(365):
            for hour in range(24):
                current_time = start_date + timedelta(days=day, hours=hour)
                time_points.append(current_time)
        
        print(f"✅ 生成时间点: {len(time_points)} 个 (全年每小时)")
        return time_points
    
    def calculate_solar_position(self, time_points):
        """计算太阳位置 (方位角和高度角)"""
        solar_positions = []
        
        for time_point in time_points:
            # 计算儒略日
            julian_day = self._julian_day(time_point)
            
            # 计算太阳赤纬角
            declination = self._solar_declination(julian_day)
            
            # 计算时角
            hour_angle = self._hour_angle(time_point, julian_day)
            
            # 计算太阳高度角
            elevation = self._solar_elevation(declination, hour_angle)
            
            # 计算太阳方位角
            azimuth = self._solar_azimuth(declination, hour_angle, elevation)
            
            # 只保留太阳高度角大于0的时间点 (白天)
            if elevation > 0:
                solar_positions.append({
                    'time': time_point,
                    'elevation': elevation,  # 高度角 (弧度)
                    'azimuth': azimuth,      # 方位角 (弧度)
                    'declination': declination,
                    'hour_angle': hour_angle
                })
        
        print(f"✅ 太阳轨迹计算完成: {len(solar_positions)} 个有效时间点")
        return solar_positions
    
    def _julian_day(self, time_point):
        """计算儒略日"""
        a = (14 - time_point.month) // 12
        y = time_point.year + 4800 - a
        m = time_point.month + 12 * a - 3
        
        jdn = time_point.day + (153 * m + 2) // 5 + 365 * y + y // 4 - y // 100 + y // 400 - 32045
        jd = jdn + (time_point.hour - 12) / 24.0
        
        return jd
    
    def _solar_declination(self, julian_day):
        """计算太阳赤纬角"""
        # 简化赤纬角计算 (Cooper方程)
        day_of_year = (julian_day - 2451545.0) % 365.25  # 年积日
        declination = math.radians(23.45) * math.sin(math.radians(360 * (284 + day_of_year) / 365))
        return declination
    
    def _hour_angle(self, time_point, julian_day):
        """计算时角"""
        # 简化时角计算 (标准公式)
        local_solar_time = time_point.hour + time_point.minute/60.0
        
        # 时角 = 15° × (当地太阳时 - 12)
        hour_angle = math.radians(15 * (local_solar_time - 12))
        return hour_angle
    
    def _solar_elevation(self, declination, hour_angle):
        """计算太阳高度角"""
        lat_rad = math.radians(self.latitude)
        elevation = math.asin(
            math.sin(lat_rad) * math.sin(declination) + 
            math.cos(lat_rad) * math.cos(declination) * math.cos(hour_angle)
        )
        return elevation
    
    def _solar_azimuth(self, declination, hour_angle, elevation):
        """计算太阳方位角"""
        lat_rad = math.radians(self.latitude)
        
        azimuth = math.atan2(
            math.sin(hour_angle),
            math.cos(hour_angle) * math.sin(lat_rad) - math.tan(declination) * math.cos(lat_rad)
        )
        
        # 转换到0-2π范围
        if azimuth < 0:
            azimuth += 2 * math.pi
            
        return azimuth
    
    def generate_candidate_positions(self, grid_size=2.5):
        """生成候选位置集合 (高密度极坐标网格) - 标准化版本"""
        candidate_positions = []
        
        # 标准化参数: 更密集的网格以确保足够的候选位置
        min_radius = self.R_keepout + 5   # 禁区外5m (更紧凑)
        max_radius = self.R_field - 5     # 场地内5m (更紧凑)
        
        radius_step = grid_size  # 2.5m径向步长 (更密集)
        base_angle_step = math.radians(2.5)  # 2.5度角度间隔 (更密集)
        
        for r in np.arange(min_radius, max_radius, radius_step):
            # 动态调整角度步长 (确保周向密度合理)
            # 内圈: 更密集; 外圈: 适度密集
            circumference = 2 * math.pi * r
            desired_spacing = grid_size * 1.2  # 期望的弧长间距
            num_points = max(8, int(circumference / desired_spacing))  # 最少8个点
            current_angle_step = 2 * math.pi / num_points
            
            for i in range(num_points):
                theta = i * current_angle_step
                x = r * math.cos(theta)
                y = r * math.sin(theta)
                
                # 检查是否在场地内且不在禁区内
                if self._is_valid_position(x, y):
                    candidate_positions.append((x, y))
        
        # 添加一些随机扰动位置以增加多样性
        random_positions = self._generate_random_positions(len(candidate_positions) // 10)
        candidate_positions.extend(random_positions)
        
        print(f"✅ 候选位置生成完成: {len(candidate_positions)} 个位置 (标准化高密度版本)")
        print(f"   径向步长: {grid_size}m, 有效半径: {min_radius:.1f}-{max_radius:.1f}m")
        return candidate_positions
    
    def _is_valid_position(self, x, y):
        """检查位置是否有效 (在场地内且不在禁区内)"""
        distance = math.sqrt(x*x + y*y)
        
        # 不在禁区内 且 在场地内
        return distance >= self.R_keepout and distance <= self.R_field
    
    def _generate_random_positions(self, num_random):
        """生成随机扰动位置以增加候选位置多样性"""
        random_positions = []
        max_attempts = num_random * 10  # 防止无限循环
        attempts = 0
        
        while len(random_positions) < num_random and attempts < max_attempts:
            # 在环形区域内随机生成位置
            r = np.random.uniform(self.R_keepout + 5, self.R_field - 5)
            theta = np.random.uniform(0, 2 * math.pi)
            
            x = r * math.cos(theta)
            y = r * math.sin(theta)
            
            if self._is_valid_position(x, y):
                # 添加小幅随机扰动 (±1m)
                x += np.random.uniform(-1, 1)
                y += np.random.uniform(-1, 1)
                
                if self._is_valid_position(x, y):
                    random_positions.append((x, y))
            
            attempts += 1
        
        return random_positions
    
    def calculate_dni_average(self, solar_positions):
        """计算年平均DNI (基于气象数据预处理)"""
        # 简化模型: 基于太阳高度角的DNI估算
        dni_values = []
        
        for pos in solar_positions:
            elevation_deg = math.degrees(pos['elevation'])
            
            # DNI与太阳高度角的经验关系 (W/m²)
            if elevation_deg > 10:
                dni = 900 * math.sin(pos['elevation'])  # 晴天模型
                dni *= 0.7  # 考虑云层等因素的年平均修正
            else:
                dni = 0
            
            dni_values.append(dni)
        
        dni_avg = np.mean(dni_values) if dni_values else 0
        print(f"✅ 年平均DNI计算完成: {dni_avg:.1f} W/m²")
        
        return dni_avg, dni_values
    
    def precompute_geometry(self, candidate_positions, tower_position=(0, 0)):
        """预计算几何关系 (定日镜-吸收塔距离等)"""
        geometry_data = {}
        
        tower_x, tower_y = tower_position
        
        for i, (x, y) in enumerate(candidate_positions):
            # 计算定日镜到吸收塔的距离
            distance = math.sqrt((x - tower_x)**2 + (y - tower_y)**2)
            
            # 计算定日镜到吸收塔的方位角
            azimuth = math.atan2(y - tower_y, x - tower_x)
            if azimuth < 0:
                azimuth += 2 * math.pi
            
            geometry_data[i] = {
                'position': (x, y),
                'distance': distance,
                'azimuth': azimuth
            }
        
        print(f"✅ 几何关系预计算完成: {len(geometry_data)} 个位置")
        return geometry_data

# 测试数据预处理模块
if __name__ == "__main__":
    print("🚀 开始数据预处理...")
    
    # 初始化预处理器
    preprocessor = DataPreprocessor()
    
    # 加载场地数据
    x_coords, y_coords = preprocessor.load_boundary_data()
    
    # 生成时间点 (测试用，只取前100个小时)
    time_points = preprocessor.generate_time_points()[:100]
    
    # 计算太阳轨迹
    solar_positions = preprocessor.calculate_solar_position(time_points)
    
    # 生成候选位置
    candidate_positions = preprocessor.generate_candidate_positions()
    
    # 计算平均DNI
    dni_avg, dni_values = preprocessor.calculate_dni_average(solar_positions)
    
    # 预计算几何关系
    geometry_data = preprocessor.precompute_geometry(candidate_positions)
    
    print("✅ 数据预处理完成！") 